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Abstract 

In this article we introduce model-independent data analysis procedures for identifying 
inelastic WIMP-nucleus scattering as well as for reconstructing the mass and the mass 
splitting of inelastic WIMPs simultaneously and separately. Our simulations show that, 
with 0(50) observed WIMP signals from one experiment, one could already distinguish 
the inelastic WIMP scattering scenarios from the elastic one. By combining two or more 
data sets with positive signals, the WIMP mass and the mass splitting could even be 
reconstructed with statistical uncertainties of less than a factor of two. 



1 Introduction 



Different astronomical observations and measurements indicate that more than 80% of all matter 
in the Universe is dark (i.e., interacts at most very weakly with electromagnetic radiation and 
ordinary matter). The dominant component of this cosmological Dark Matter (DM) must be 
due to some yet to be discovered, non-baryonic particles. Weakly Interacting Massive Particles 
(WIMPs) x arising in several extensions of the Standard Model of electroweak interactions are 
one of the leading DM candidates. Currently, the most promising method to detect different 
WIMP candidates is the direct detection of the recoil energy deposited in a low-background 
underground detector by scattering of ambient WIMPs off target nuclei (for reviews, see Refs. [1, 
2, 3, 4, 5, 6]). ^ _ ^ 

The first positive signal of WIMPs was reported by the DAMA Collaboration with obser- 
vation of the annual modulation of the event rate of possible DM-target interaction by the 
DAMA/Nal detector [7, 8] (for the latest updated results, see Ref. [9]). The CoGeNT Collabo- 
ration announced also their result with a modulated component of unknown origin [10, 11]. The 
CRESST Collaboration published in 2011 their newest result with more observed events than 
expected backgrounds in the CRESST-II experiment [12]. However, other direct DM detection 
experiments can so far mostly observe only very few candidate events and a large part (~ 40% 
to 60%) of these events would be unrejected backgrounds [13, 14, 15, 16, 17]. 

Many theoretical scenarios have been proposed to find an explanation for reconciling all 
these results from different experiments with various target nuclei. Among them, one of the 
mostly discussed and experimentally constrained frameworks is inelastic Dark Matter (iDM) 
models [18, 19]. The basic assumption of iDM models is that the incident WIMPs scatter 
inelastically off target nuclei and then transit to a slightly heavier (exciting) state with a tiny 
mass (energy) splitting 5, while the elastic scattering channel is suppressed or forbidden. Due to 
the kinematics of inelastic scattering, the targets with heavier nuclei have greater sensitivities 
than those with lighter nuclei. This suggests that the signals reported by DAMA which utilises 
relatively heavy iodine nuclei may not be detected by other experiments with light target nuclei 
such as germanium and silicon. Comparisons between exclusion limits on the cross section versus 
mass (splitting) (a — m x {5)) planes in the iDM scenarios by using various experimental data 
have been done [20, 21]. In Ref. [21] the authors considered also the case of low-mass WIMPs 
(m x ~ 5 GeV) with a tiny mass splitting (5 ~ 10 keV). Meanwhile, a method for reconstructing 
the iDM parameters, most importantly the WIMP mass m x and the mass splitting 5, based on 
likelihood analysis has been suggested [22]. 

In the recent years, several experimental collaborations have (re) analyzed their (null) obser- 
vation results and found severe constraints on the parameter space of iDM, i.e., on the mass 
splitting 5 versus the WIMP mass m x plane. From these analyses, exclusion limits for mass 
splitting up to 250 keV in the mass range between 40 GeV and 1 TeV with different detector 
materials have been given [13, 23, 24, 25, 26, 17]. However, in all these works, data analyses 
have been done basically by estimations of the (differential) event rate, which is strongly halo- 
model dependent. In addition, types and (relative) strengths of different DM-nucleus (quark) 
interactions are important factors in direct and indirect DM detection experiments [1, 2]. 

Hence, in this paper, as complementarity and extension of our earlier work on develop- 
ing methods for reconstructing WIMP properties in elastic scattering framework as model- 
independently as possible [27, 28, 29, 30], we introduce new model-independent approaches for 
identifying inelastic WIMP-nucleus scattering scenarios as well as for reconstructing the most 
important properties of inelastic WIMPs: the mass m x and the mass splitting S simultaneously 
and separately. Our method is based on an estimation of a characteristic energy Q VthIC corre- 



sponding to a threshold (minimal required) (one-dimensional) velocity of incident (inelastic) 
WIMPs, fthre- 1 This can be done by determining the maximum of the integral over the one- 
dimensional velocity distribution function of incident WIMPs. Once this characteristic energy 
can be solved by using data sets of positive (inelastic) WIMP signals with different target nuclei, 
one can then determine the (degenerate) WIMP mass and the tiny mass splitting straightfor- 
wardly. 

The remainder of this article is organized as follows. In Sec. 2 We develop the formalism 
of the methods for reconstructing the one-dimensional WIMP velocity distribution function, 
determining the WIMP mass and the mass splitting as well as estimating the characteristic 
energy and in turn identifying the inelastic WIMP-nucleus scattering scenarios. In Sec. 3 we 
demonstrate the ability and shortcomings of our model-independent methods by presenting nu- 
merical results based on Monte-Carlo simulations. The possibility of distinguishing the inelastic 
WIMP scenarios from the elastic one will be particularly discussed. We conclude in Sec. 4. Some 
technical details for our analysis will be given in the appendices. 



2 Formalism 



In this section we develop the formulae needed in our model-independent reconstructions of 
different properties of halo WIMPs in iDM scenarios, both of an approximated analytic method 
and an iterative numerical procedure will be considered. 

We start with the basic expression for the differential event rate for (elastic) WIMP-nucleus 
scattering given by [1]: 



fi(v) 



dv . 



(1) 



Here R is the direct detection event rate, i.e., the number of events per unit time and unit mass 
of detector material, Q is the energy deposited in the detector, F(Q) is the (elastic) nuclear form 
factor, /i(i>) is the one-dimensional velocity distribution function of the WIMPs impinging on 
the detector, v is the absolute value of the WIMP velocity in the laboratory frame. The constant 
coefficient A is defined as 



A 



2m x m? N 



(2) 



where p is the WIMP density near the Earth and a is the total cross section ignoring the form 
factor suppression. The reduced mass m T N is defined by 



m riN 



(3) 



where m x is the WIMP mass and that of the target nucleus. Finally, t> min is the minimal 
incoming velocity of incident WIMPs that can deposit the energy Q in the detector: 



y/2m N Q 



m r , N 



Q + 5 



ctJQ + 



a 5 



(4) 



1 Remind that, in conventional elastic scattering framework, all incident WIMPs with non-zero velocity could 



scatter off target nuclei and deposit recoil energies Q. This means that Q v 



for clastic WIMPs and 



Qvthre > for inelastic WIMPs. In fact, later we will show that the characteristic energy Q Vthre is proportional 
to the mass splitting 5. 



with the transformation constants 



and 

t) mEH is the maximal WIMP velocity in the Earth's reference frame, which is related to the escape 
velocity from our Galaxy at the position of the Solar system, t> esc > 600 km/s. 



2.1 Reconstruction of the velocity distribution fi(v) 

We first define [27] 

dF^v) _ fjv) 
dv v ' 

i.e., Fi(v) is the primitive of fi(v)/v. Eq. (1) can then be rewritten as 2 

' ' dv = F ± (v = f max ->■ oo) - 



(7) 



(8) 

Since WIMPs in today's Universe move quite slow, /i(i>) must vanish as v approaches infinity, 
./i(u->oo)->0. (9) 
Hence 



/d-R\ /-fmax->00 


r/i(«)i 


Y (iQ / ■/ f min 





-> 0. 



(10) 



This means that Fi(t> — >■ oo) approaches a finite value. Differentiating both sides of Eq. (8) with 
respect to i> m i n and using Eq. (4), we can find 



dFi(t> min ) 
dv m i n 



d 



1 



v4 I dv m i n 
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as 



-2Q 



= Q(fmin) , 

d 



dg 



'dR s 



F 2 (Q) \dQ / 



(11) 



Q=Q(Vn 



Here, from the definition (4) of f m in, we can firstly find that 

d^rnin 1 ( Ct5 \ 

aJQ - 



dQ 2Q 



VQ ' 



;i2) 



and an analytic expression of Q(v m m) can be solved from the definition (4) directly as: 

Q (^min 



W min ~ 2aa S ± fmin\/^min _ ^ aa S 



2a 2 



(13) 



2 For simplicity, we set at first here the maximal cut-off of the one-dimensional WIMP velocity distribution 
as infinity. Later we will discuss correction of the formulae given here in the practical use with real (generated) 
data events. 



Note that, corresponding to one specific value of i> m i n , there are two possible values of Q(v min ), 
except of 



Qv thla = Q (Vmin = fthrc = 2y/aa~s~) = — = f ^ 

a \m x + m N/ 



5. 



(14) 



Here i> t h re is the threshold (minimal required) velocity of incident inelastic WIMPs, which can 
produce recoil energy at all. 

Since the expression (11) of dFi(v m i n ) / dv m i n holds for arbitrary t> m in, we can write down the 
following result directly: 



fi(v) dF 1 (v) 



dv 



as 
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'dR\ 



dQ [F 2 (Q) \dQ) 



(15) 



Q=Q(v) 



with 



Q(v) 



v 2 — 2aas ± v \/v 2 — Aacts 
2^ 



(16) 



Although the right-hand side of this expression depends on the as yet unknown constant A, 
fi{v) is the normalized velocity distribution, i.e., it is defined to satisfy 



fi(v) dv=l. 



(17) 



Therefore, the normalized one-dimensional velocity distribution function of inelastic WIMPs 
can be given by 



f 1 ( v )=Ml\aJQ + 



as_ 



aJQ — 



as 



with the normalization constant M: 



M = 2 



as 



dR' 



[F 2 (Q) \dQj\ 
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F 2 (Q) UQ / 



,(18) 



=<?(«) 



(19) 



Remind that, for the case of elastic WIMP-nucleus scattering, as = 0, the expressions (18) and 
(19) can then be reduced to the simple, analytic forms given in Eqs. (12) and (13) of Ref. [27]. 



2.2 Determinations of the WIMP mass m x and the mass splitting 8 

The expression for reconstructing the one-dimensional velocity distribution of inelastic WIMPs 
given in Eqs. (18) and (19) can unfortunately not be used directly, since at first there are two 
unknowns, i.e., the WIMP mass m x (involved in a) and the mass splitting 5 (involved in as)- 
Moreover, the typical peaky shape of the recoil spectrum of inelastic WIMPs (see e.g., Figs. 1 to 
4) makes its reconstruction more complicated and therefore a similar development of procedures 
introduced in Refs. [27, 28] is basically impossible. 

Hence, in this subsection we introduce a new approach for determining the WIMP mass and 
the mass splitting based on the estimation of the characteristic energy Q Vthia , which requires the 
reconstruction of the peaky inelastic WIMP recoil spectrum. 



2.2.1 Determinations of m x and S 

From the definition (4) of i> m i n , one have 



dvr. 



dQ 



Q=Q 
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= 0. 



(20) 



Then the characteristic energy corresponding to the minimal value of i> m i n , fthrc, can easily be 
solved as 



as 



5, 



(14) 



a \m x + m^ / 

which is proportional to the mass splitting 5 and the proportionality constant is simply a function 
of the WIMP mass m x . Hence, by combining two experimental data sets with different target 
nuclei, X and Y, one can derive analytic expressions for determining m x and 5 as functions of 
the characteristic energy Q Vthiet x and Q Vthie , Y - 

Qv thIC ,Ym Y - Q VthTe ,xm x 



™x = 
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Q 



and 



5 = 



Qv thTe ,xQv thTC ,Y (m Y - m x ) 



(21) 



(22) 



Qv thlc ,Y m Y ~ Qv thic ,xm x 

Then by using the standard Gaussian error propagation, the statistical uncertainties of m x and 
5 can be given as 
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and 



a (5) 
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(24) 



2.2.2 Ansatz for reconstructing the inelastic-scattering recoil spectrum 

Remind that, for the use of Eqs. (21) and (22), one need to estimate the characteristic energy 
Qvthre corresponding to the threshold (minimal required) velocity of incident inelastic WIMPs 
Vthie, which could produce recoil energy at all. This means that v = v t hre is the lowest bound of 
the velocity of inelastic WIMPs, which could contribute to the integral over the one-dimensional 
velocity distribution function fi(v) on the right-hand side of Eq. (1). This means in turn that 
the integral over fi(v), or, equivalently, a "reduced" differential event rate (i.e., the differential 
event rate divided by the squared nuclear form factor): 



'dR s 



F 2 (Q) UQ / 
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fi(v) 



dv 



(25) 



should be maximal once Q = Q Vthie - 

Now, for the simplest isothermal spherical halo model, the normalized one-dimensional ve- 
locity distribution function can be expressed as [27]: 



v 2 \ ,J ■ ._> 



h ^ = T^V~ vK ' (26) 

where vq ~ 220 km/s is the Sun's orbital speed around the Galactic center. On the other hand, 
by taking into account the orbital motion of the Solar system around the Galaxy, the more 
realistic, and thus more frequently used shifted Maxwellian velocity distribution has been given 
by [27]: 

1 / v 

where v e is the time-dependent Earth's velocity in the Galactic frame: 



/l,Bh(u) 



e -(v-V e ) 2 /v$ _ e -(f+t>c) 2 A>0 



(27) 



v e (t) = v 



1.05 + 0.07 cos 



'2n(t-t p ) 

. 1 y r 



(28) 



where t p ~ June 2nd is the date on which the Earth's velocity relative to the WIMP halo is 
maximal. 

Substituting these two expressions into Eq. (1) and using Eq. (4), one can easily obtain that 3 

(^) oc e A-VQ +as /VQ) 2 /^ K e -{^Q +a y Q )/ v i ? (29a) 

F (Q) \dQJ Qau 



and 



' f dB ?\ oc erf ("v / Q+^/VQ+^) _ eii ^VQ+^/VQ-v^ _ (29b) 



F 2 (Q) UQ 



in, sh 



Then, similar to the use of the exponential approximation for reconstructing the recoil spectrum 
of elastic WIMP-nucleus scattering [27], in order to approximate the measured recoil spectrum 
and take into account the extra contribution predicted in inelastic WIMP scattering scenarios, 
we introduce empirically here 

(|). (30) 

\ ^ / m, cxpt 



where 

r o = fQm „ — . , ( 31 ) 



£ I r m - e -*Q-*/Q dQ 



with the total WIMP signal events in our data set JV to t and the experimental exposure £; k and 
k' are two fitting parameters which we want to estimate by using the measured recoil energies 
directly. By using the ansatz (30), the position of the peak of the inelastic recoil spectrum can 
be solved easily and then estimated by (mathematical details are given in Appendix A. 2) 

[Q- l '% 

)inf 



~ /g-3/2\ mf ' ( 32 ) 



3 As in Sec. 2.1, we assume here that the cut-off on fi(v), f max , as well as the experimental maximal cut-off 
energy Q max are large enough and the integral in the higher velocity/energy range can be neglected. 



with the statistical uncertainty given by 
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(33) 



where the A— th momentum of the recoil energy spctrum can be estimated by 
Q x (dR/dQ) dQ i 

V / in, cxpt 1 



Q x ) 



inf 
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tot a 



(34) 



Note that we assumed here (unrealistically) that the minimal experimental cut-off energy Qmin 
is negligible and the maximal one Q max is infinity (a numerical correction will be discussed in the 
next subsection). Finally, two fitting parameters in the ansatz (30) can be estimated separately 

as 
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and 
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(35b) 



By using the approximation (30), the expression (18) for reconstructing the one-dimensional 
WIMP velocity distribution can be rewritten as 



fi(v)=tfl[aJQ + 
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(36) 



Q=Q(v) 



with the normalization constant M given in Eq. (19). Note that, for reconstructing f\(v) by 
using Eqs. (36) and (19), the constant r in Eq. (31) can be cancelled out, since this appears in 
both of the expression (18) or (36) and the expression (19). 

On the other hand, as discussed at the beginning of this subsection, i> thre is the lowest bound 
of the integral in Eqs. (1) or (25), which gives a maximal value of the (reduced) event rate and 
thus have to satisfy the following condition: 
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F 2 (Q) \dQ / 

By using the ansatz (30) for reconstructing the inelastic scattering spectrum, one can find that 



(38) 



Q Vth can then be solved numerically and the statistical uncertainty on Q Vth can be estimated 

by " 



a{Q 



Vthr. 





E 

A,p=-3 



'dQ 



dk 



dk 



d{Q x+1 ' 2 ) : 



+ 



'dQ 



- '-'thr. 



dk' 



dk' 



x 



'9g 



i>thr. 



dk 



'dQ 



xcov((Q a+1 / 2 ),(Q" +1 / 2 



1/2 



d(Q 



p+l/2 



k! d 
Q* + dQ 



1 



'dF' 



x 



o 

E 

A,p=-3 



F(Q) \dQ) 
dk 



dk' 



\ d (Q x+1/2 ), 



Q 2 



Vthr 



d(Q 



A+l/2 



X 



dk 



dk' 



d(Q 



p+l/2 



Q 



Vthr 



X cov 



((Q A+1/2 ),(Q' +1/2 



d{Q p+l ' 2 

s 1/2 



Here we have 
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(40) 



and dk^ / d (Q x+1/2 ) for A = -3, -2, - 1, are given in Appendix C. Note that, for the 
analytic estimates of k and k' given by Eqs. (35a) and (35b), one has 



dh 



dk' 



= 0. 



(41) 



2.3 A numerical correction 



By using Eqs. (34), (35a) and (35b), one has to assume that the (experimental) minimal cut-off 
energy should be negligibly small (Q m m — 0) and the maximal one large enough (Q max ^ a few 
hundred keV). This should be an achievable requirement for new generation detectors. How- 
ever, due to the local escape velocity of halo WIMPs, v esc > 600 km/s, or, equivalently, the 
maximal cut-off on the one-dimensional WIMP velocity distribution, f max , there are not only a 
kinematical maximal but also a minimal cut-off energies, with which the incident WIMPs could 
produce recoil energies in our detector (see e.g., Figs. 4). From Eq. (16), one can get 
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(42b) 



As shown later in Sec. 3, this could be a serious issue once the mass splitting 5 > 30 keV, 
especially for light target nuclei, i.e., 28 Si and 40 Ar (see Figs. 4). 

Moreover, as we will show in Sec. 3.3, for larger WIMP masses (m x > 200 GeV), the recon- 
structed WIMP mass could be strongly underestimated: the larger the mass splitting, the worse 
the mass reconstruction. Hence, we introduce in this subsection an iterative method for cor- 
recting the estimation of the fitting parameters, k and k', in the hope that this correction could 
reduce the systematic deviation of the solved Q VthIC and thereby alleviate the underestimation 
of the reconstructed WIMP mass. 

We start from Eq. (34). Define 

^(^- ^M^Z^M^^ , (43 ) 

" tc/m i n 



where (x J (k, k'; x) are given in Eqs. (A17) to (Al8c) in Appendix A. 2. In stead of \Q , . 

defined on the left-hand side of Eq. (34), (<? A ) (k, k') defined above should be more suitable to 
be estimated by the average value of the summation over all measured recoil energies in the Ath 
power given on the right-hand side of Eq. (34). Thus one has 



Qr, 



(44) 



By setting fcW* = k^ a estimated by Eqs. (35a) and (35b), one can solve fcW numerically, denoted 
by k^l m , by using the expressions of (x x ) (k,k';x) , A = — 1/2 and —3/2, simultaneously 4 . 

» ' Qmin 

Note that expressions (36) and (38) for reconstructing the one-dimensional WIMP velocity 
distribution and solving the characteristic energy Qt, thrc can be used by substituting both the 
analytic and numerical estimations of the fitting parameters k and k', whereas for using Eq. (39) 
to estimate the statistical uncertainty on Q VthIC solved with the numerically estimated k num and 
k' nura , the summation runs only for X, p = —2 and —1. 



3 Numerical results 

In this section, we present numerical results of the reconstruction of the recoil spectrum, the 
identification of the positivity of Q Vtbrc (as the check of the iDM scenarios) as well as the recon- 
struction of the WIMP mass m x and the mass splitting 5 based on Monte-Carlo simulations. 
The special case of zero mass splitting 5 = 0, i.e., the case of elastic WIMP scattering, will 
be particularly discussed at the end of this section as a demonstration of the usefulness of our 
model-independent approach for distinguishing the inelastic WIMP scenarios from the elastic 
one. 

For generating WIMP-induced signals in our simulations, we use the shifted Maxwellian ve- 
locity distribution given in Eq. (27) with the Sun's Galactic orbital velocity t>o — 220 km/s; the 
time dependence of the Earth's velocity in the Galactic frame has been ignored, i.e. v e = 1.05 vo 
has been used. Moreover, the maximal cut-off of the one-dimensional WIMP velocity distribu- 
tion has been set as v max = 700 km/s. Meanwhile, the WIMP-nucleon cross section has been 



Qmax 

'The detailed description about solving fc nU m and k' num from (a; -1 / 2 ) (k, k'; x) and (a; -3 / 2 ) (k, k'; x) 



simultaneously will be given in Appendix B. 



mm 



assumed to be only spin-independent (SI), a^ p = 10 
for the elastic nuclear form factor: 



pb, and the commonly used analytic form 



(Q) 



3ji(g-Ri 
qRi 



(45) 



has been adopted. Here Q is the recoil energy transferred from the incident WIMP to the target 
nucleus, ji(x) is a spherical Bessel function, q = \/2m^Q is the transferred 3-momentum, for 
the effective nuclear radius we use R-y = y/R\ - 5s 2 with R A ~ 1.2 A 1/3 fm and a nuclear skin 
thickness s ~ 1 fm. In addition, the experimental threshold energies have been assumed to 
be negligible (Q m in = 0) and the maximal cut-off energies are set as Q ma , x = 150 keV for all 
target nuclei. 5,000 experiments with 50 total events on average in one experiment have been 
simulated. 



3.1 Reconstructing the recoil spectrum 

In this subsection, we consider at first the reconstruction of the recoil spectrum, which is ap- 
proximated by Eq. (30) with the fitting parameters k and k', as well as the estimation of the 
characteristic energy Q VthIC - 

The top-left frame of Figs. 1 shows the measured recoil energy spectrum (dotted magenta 
histogram) for a 76 Ge target. The dash-double-dotted cyan curve is the inelastic WIMP-nucleus 
scattering spectrum used for generating recoil events. The input WIMP mass and mass splitting 
are set as m x = 100 GeV and 5 = 25 keV. We show here also the reconstructed recoil spectra: 
while the dashed blue curve is the 2-parameter exponential spectrum (30) with the parameters 
k and k! estimated analytically by Eqs. (35a) and (35b) 5 , the solid red and dash-dotted green 
curves are given with the parameters k and k' estimated by the numerical method introduced 
in Sec. 2.3; the former and later show the results obtained from the first and final rounds of the 
iterative process 6 . 

It can be found that, firstly, the recoil spectra reconstructed numerically could really be the 
better approximations of the original theoretical spectrum, at least in the sense that the peaks of 
these spectra coincide very well. Secondly, and not really as expected, although the numerically 
reconstructed recoil spectra could be the better results than the analytically reconstructed one, 
the iterative procedure could not improve the correction further: only in a small part of the 
simulated experiments (shown here and later with different initial setup) the results coming 

5 Remind that, for the analytically reconstructed spectrum shown here, the constant r in Eq. (30) has been 
estimated by Eq. (31), in which the lower and upper bounds of the integral in the denominator have to be set as 
Qmin = and Q m ax = oo. Then Eq. (31) can be rewritten as 

(46) 

where K\{x) is the modified Bessel function of the second kind of order 1. Detailed derivation can be found in 
Appendix A.l. 

6 Note that, as we will discuss later, the results given by our iterative numerical procedure might diverge 
(deviate larger and larger from the theoretical values) round by round. Thus, once the results given by the nth 
run is too far away from the (n — l)th run, our program will take the results of (n — l)th run as the final results. 
This means that in some of the 5,000 simulated experiments the results of the final round are just those of the 
first or the second rounds. 

Note also that, because of the same divergency problem and of that the results given by the later round might 
not be much better than those of the first round, for the reconstructions of Q Vthro as well as the WIMP mass and 
the mass splitting shown later, we use only the results (Lm and k' num ) given by the first round. 
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Figure 1: The measured recoil energy spectrum (dotted magenta histogram) for a 76 Ge target 
(top-left) as well as the distributions of the reconstructed fitting parameters k (top-right), k' 
(bottom-left) and the characteristic energy Q VthTC (bottom-right). The dash-double-dotted cyan 
curve is the inelastic WIMP-nucleus scattering spectrum used for generating recoil events. The 
input WIMP mass and mass splitting are set as m x = 100 GeV and 5 = 25 keV. While the 
dashed blue curve is the 2-parameter exponential spectrum (30) with the parameters k and k' 
estimated analytically by Eqs. (35a) and (35b) (top-left) or the distributions of the analytically 
estimated k, k' and Q Vthrc , the solid red and dash-dotted green curves are given with or for the 
parameters k and k' estimated numerically; the former and later show the results obtained from 
the first and final rounds of the iterative process. Meanwhile, the cyan (magenta, light green) 
vertical lines indicate the median values of the simulated results (corresponding to the blue, red 
and green histograms), whereas the horizontal thick (thin) bars show the 1 (2) a ranges of the 
results. The green vertical line given additionally (bottom-right) indicates the theoretical value 
of Q VthIC estimated by Eq. (14), Q Vthrc ,th- See the text for further details. 

from the second (and also the final) rounds could be clearly better than those coming from 
the first round; in most part of simulations there would no (significant) difference between the 
results given by different rounds, or those from the later round could even be worse... 

Moreover, in the top-right and bottom-left frames of Figs. 1 we show the distributions of 
the reconstructed fitting parameters k and k', respectively: while the dashed blue curves show 
the distributions of k and k' estimated analytically by Eqs. (35a) and (35b), the solid red and 
dash-dotted green curves are those of the numerically estimated k and k', the former and later 
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Figure 2: As in Figs. 1, except that 28 Si (top), 40 Ar (middle) and 136 Xe (bottom) have been used 
as target nuclei. 



show the results given by the first and final rounds of the iterative process. Meanwhile, the 
cyan (magenta, light green) vertical lines indicate the median values of the simulated results 
(corresponding to the blue, red and green histograms), whereas the horizontal thick (thin) bars 
show the 1 (2) a ranges of the results 7 . 

The top-right frame of Figs. 1 shows that the distributions of the reconstructed fitting 
parameter k by both of the analytic and numerical methods are indeed basically Gaussian (with 
small tails in the high-fc range). Meanwhile, in the bottom-left frame the distributions of 
the reconstructed k! by both methods would be Gaussian distributions with a k' = cut-off. 
Remind that, k' is basically proportional to the square of the mass splitting, 5 2 , (see Eq. (29a) 
and then Eq. (6)). The observation of the positively reconstructed k' indicates that our model- 
independent reconstruction proposed here should be useful and sensitive for identifying the 
inelastic WIMP-nucleus scattering scenarios. 

Moreover, these two plots show also that there are neither significant differences between 
the median values of k and k' (and thus the reconstructed recoil spectra) estimated in the first 
and final rounds of the numerical procedure, nor those between the distributions of them in 

7 Here the 1 (2) a ranges mean that 68.27% (95.45%) of the reconstructed values in the simulated experiments 
are in these ranges (central interval). 
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Figure 3: As in Figs. 1 and 2, except that a smaller mass splitting 5 = 10 keV has been used. 

the simulated experiments. In contrast, the differences between the numerical results and the 
analytic ones are clear. Hence, later we show only results given by the first round of the iterative 
process as the numerical reconstruction. 

Finally, in the bottom-right frame of Figs. 1, we show the distributions of the characteristic 
energy Q VthIB solved by Eq. (38) with k and k! estimated analytically (dashed blue) and nu- 
merically (solid red), respectively. The green vertical line given additionally here indicates the 
theoretical value of Q Vthia estimated by Eq. (14), Q Vthie ,th- This plot shows that, firstly, both of 
the distributions of Q VthIC solved analytically and numerically are still basically Gaussian. Sec- 
ondly, the median value of the numerically solved Q Vthre is indeed much closer to the theoretical 
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estimate with a smaller 1 (2) o upper and lower statistical uncertainties (and however a little bit 
longer tail in the high-energy range). 

Considering the currently running and built/planned next-generation detectors, in Figs. 2 
we present the results with several different target nuclei: 28 Si (top), 40 Ar (middle) and 136 Xe 
(bottom). Basically, all observations discussed above hold, except that the median value of the 
analytically solved Q VthIC with the 136 Xe target is now closer to the theoretical estimate. 

In Figs. 3 we consider a smaller mass splitting 5 = 10 keV. For this extreme case 8 , firstly, 



8 The special case of 5 = (elastic WIMP-nucleus scattering) will be discussed particularly at the end of this 



the tails of the distributions of k, k' and Q VtbTC solved analytically and numerically found in 
Figs. 1 and 2 are reduced significantly. This indicates that, for smaller mass splittings, the 
statistical fluctuations of k, k' and Q Vthie given by both of the analytic and numerical methods 
would in principle be smaller. However, while there is still no significant difference between the 
numerically estimated k' from different rounds of the iterative procedure (the third column), the 
estimations of the fitting parameter k (the second column) become clearly worse (and worse). 
This causes in turn larger deviations of the reconstructed recoil spectrum (see the first column). 

Moreover, in all four plots of the distribution of the solved Q Vthie (forth column), an excess 
around Q Vthia ~ can be seen clearly. This is caused by our setup of a lower bound used 
in solving Q Vthie and means thus that there is a small (but non-zero) possibility of obtaining 
non-physically negative Q Vthia , once the mass splitting is pretty small. 

On the other hand, we would like to remind here that, as shown in the first and forth columns 
of Figs. 3, since Q VthIC is proportional to the mass splitting S (for a fixed WIMP mass m x , see 
Eq. (14)), once the mass splitting is quite small (5 < 10 keV), the position of the peak of the 
recoil spectrum, Q pk (not, but related to Q VthIC ), would also be pretty small, probably smaller 
than the experimental minimal (software or even hardware) cut-off energy 9 . Hence, a reduction 
of the threshold energy to be small enough/negligible would be necessary. Nevertheless, once this 
requirement can be achieved, our simulations shown here and in the next subsections indicate 
that a model-independent identification of inelastic WIMPs would in principle be possible, even 
for a small mass splitting of 10 or 5 keV. 

In contrast, in Figs. 4 a larger mass splitting 5 = 50 keV has been used. For this (relatively) 
high-5 case, the effect of the kinematic minimal and maximal cut-off energies (given in Eqs. (42a) 
and (42b)) becomes serious, especially for lighter target nuclei, e.g., 28 Si (second row). In this case 
the analytic results seem to be more reliable; the tails of the distributions of the reconstructed k, 
k' and Q VthIC offered by the numerical iterative procedure in high- (and even low-)value ranges 
become much longer and the divergency of the results mentioned earlier becomes now more 
problematic. Only numerical results offered from experiments with heavy target nuclei, e.g., 
76 Ge and 136 Xe, could be used as auxiliary. Note here that the small bump of the distribution 
of Qv thrc reconstructed with the 136 Xe nucleus (forth frame in the bottom row) between 30 and 
35 keV is caused by the (first) zero pont of the used nuclear form factor given in Eq. (45). 



3.2 Identifying the positivity of Q Vthre 

In this subsection, we study the positivity of the reconstructed characteristic energy Q VthIC (as the 
most important criterion of the identification of inelastic WIMP-nucleus scattering) in details. 

In the left frame of Figs. 5 we show the confidence levels (CLs) of the positivity of the ana- 
lytically reconstructed Q Vthie (reconstructed Q Vthie in units of the estimated la lower statistical 
uncertainties 0~\ o (Q Vthie ))' > 5er CL (blue filled circles), > 4<r CL (red filled squares), > 3a CL 
(green filled diamonds), > 2a CL (magenta filled up-triangles) , > la CL (cyan filled pentagons), 
< la CL (orange up-half-filled circles). Here we use 76 Ge as the target nucleus and check 21 
different input WIMP masses between 5 GeV and 1 TeV with 21 different input mass splittings 
between and 200 keV. Note that the empty areas on the upper part of these plots indicate 
that the ability for reconstructing the fitting parameters k and k' and then for solving Q VthIC 
would be limited by the kinematic minimal and maximal cut-off energies <3(min,max),kin given in 

section. 

9 For very light target nuclei, e.g., 19 F, Q Vthte is a little bit higher. Thus currently running and the next- 
generation detectors with such materials would be more suitable for identifying low-5 inelastic WIMPs [31, 32, 
33, 34, 35]. 




Fi gure 5: Left: confidence levels (CLs) of the positivity of the analytically reconstructed Qv thIC 
(reconstructed Q VthIO in units of the estimated la lower statistical uncertainties cr lo (Q Vthie ))' > 5cr 
CL (blue filled circles), > 4cr CL (red filled squares), > 3a CL (green filled diamonds), > 2a CL 
(magenta filled up-triangles), > la CL (cyan filled pentagons), < la CL (orange up-half-filled 
circles). Right: deviations of the analytically reconstructed Q Vthia from the theoretical values 
(A = Qv thIe ,th — Qv thvc ,rec) in units of the estimated la lower statistical uncertainties a\ Q (Qv thie )' 
2a < A (blue filled circles), la < A < 2a (red filled squares), < A < la (green filled dia- 
monds), — la < A < (magenta filled up-triangles), —2a < A < — la (cyan filled pentagons), 
A < —2a (orange up-half-filled circles). Here we use 76 Ge as the target nucleus and check 21 
different input WIMP masses between 5 GeV and 1 TeV with 21 different input mass splittings 
between and 200 keV. Other parameters are as in Figs. 1. See the text for further details. 



Eqs. (42a) and (42b), since either the WIMP mass is too small or the mass splitting is too large. 
It can be found here that, for WIMP masses < 150 GeV (with 76 Ge as the target nucleus), the 
analytically reconstructed Q Vthre should in principle be at least 5a CL apart from zero; for larger 
masses 200 GeV < m x < 1 TeV, Q VthIC should still be 3a — 4a CL apart from zero. 

However, as shown in the previous subsection, the (analytically) reconstructed Q Vthia would 
be either overestimated (for smaller 5) or underestimated (for larger 5). Hence, in the right frame 
of Figs. 5 we also check the deviations of the analytically reconstructed Q VthIC from the theoret- 




Figure 6: As in Figs. 5, except that Q Vthre have been reconstructed numerically. 
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Figure 7: As in Figs. 5 (analytically estimated Q Vthr J, except that 28 Si (top), 40 Ar (middle) and 
136 Xe (bottom) have been used as target nuclei. 

ical values (A = Q Vthrct th ~ Qv thIC ,mc) in units of the estimated la lower statistical uncertainties 
°~io {QvthrJ'- 2cr < A (blue filled circles), la < A < 2a (red filled squares), < A < la (green 



filled diamonds), — la < A < (magenta filled up-triangles) , — 2a < A < — la (cyan filled 
pentagons), A < —2a (orange up-half-filled circles). It can be found that, firstly, for WIMP 
masses < 50 GeV (with 76 Ge as the target nucleus), the analytically reconstructed Q Vthrc would 
be overestimated (A < 0), once the mass splitting 5 < 40 keV; this upper bound could be 
reduced to ~ 30 keV for heavier WIMP masses (m x > 300 GeV). Secondly, for cases of the 
non-zero mass splitting (5 > 0), the deviations would in principle be maximal 2a or less than 
la overestimated (—2a\ (Q Vthre ) < A). 

Now, by combining results showing in two frames of Figs. 5, one could conclude that, although 
the analytically reconstructed Q VthIC could be la — 2a overestimated, with around 50 total events 
one could in principle already identify the positivity of Q Vthrc with a high (3a to 5a) confidence 
level, except of the 5 = (elastic WIMP scattering) case. This observation can in turn be used 
for distinguishing the inelastic WIMP-nucleus scattering scenarios from the elastic one. 

As comparison, in Figs. 6, we show the confidence levels of the positivity as well as the 
deviations of the numerically reconstructed Q VtbTC with 76 Ge as the target nucleus. As shown in 
Figs. 1, 3 and 4, the statistical uncertainty on the numerically reconstructed Q VtbIC as well as 
Qv tbTC itself would in principle be smaller than the (uncertainty on the) analytically reconstructed 
one. Hence, firstly, for WIMP masses < 300 GeV, the numerically reconstructed Q Vthre could 
be 5a CL apart from zero; for larger masses 300 GeV < m x < 1 TeV, Qv thrc could still be 
4o" CL apart from zero. Meanwhile, in the right frame it can be seen clearly that the boundary 
line between over- and underestimations of Q VthIC is reduced to ~ 30 keV (for lighter WIMPs 
m x < 100 GeV) to ~ 20 keV (for heavier WIMPs m x > 100 GeV). 

Moreover, in Figs. 7, we check the confidence levels of the positivity as well as the deviations 
of the analytically 10 reconstructed Q Vthrc with 28 Si (top), 40 Ar (middle) and 136 Xe (bottom) as 
detector materials. Although results offered by using light target nuclei, e.g., Si and Ar, would 
(almost) always be overestimated, an at least 3a difference between reconstructed Q VtbTC and zero 
could still be observed. Meanwhile, by using heavier nuclei, e.g., Xe, one could not only test a 
wilder area on the 5 — m x plane but also give results with a higher confidence level. 

3.3 Determining m x 

In this subsection, we present the simulation results of the reconstruction of one of the key 
properties of inelastic WIMPs, i.e., the WIMP mass rn x . n 

In the left frame of Figs. 8, we show the reconstructed WIMP mass m x estimated by Eq. (21) 
and the lower and upper bounds of the la statistical uncertainty estimated by Eq. (23) as 
functions of the input WIMP mass for the case of 5 = 25 keV by using 28 Si and 76 Ge as two 
target nuclei 12 with 50 total events on average in each data set. The dashed blue curves indicate 

10 Remind that, as discussed in the previous subsection, the analytically reconstructed fitting parameters k and 
k' as well as the further solved Qt, thro would in principle be more reliable for higher mass splittings. 

11 Remind that in all simulations shown in this and the next subsection, we always reconstruct the parameters 
m x and 6 simultaneously in each simulated experiment. This means that neither m x nor 5 has been fixed (as the 
input value) in our data analysis procedure. 

12 From our results shown in Figs. 5 to 7, it seems that the WIMP mass (and the mass splitting) would be 
better reconstructed by combining two heavy target nuclei. However, our simulations show that, in practice, a 
combination of one light and one heavy target nucleus should be more suitable for reconstructing the WIMP 
mass. This can be understood as follows. The estimator given in Eq. (21) (Eq. (22)) of the (statistical uncertainty 
on the) reconstructed WIMP mass is inversely proportional to the (squared) difference between the characteristic 
energies of the recoil spectrum of two target nuclei (multiplied by the atomic mass of each nucleus). Hence, since 
the difference between the characteristic energies of the recoil spectrum of two heavy target nuclei is (much) 
smaller than that of one light and one heavy target nucleus, the statistical fluctuation caused by the use of 
only a few tens events (from one data set) would affect strongly our estimation of the (median values of the) 



Figure 8: Left: reconstructed WIMP mass and the lower and upper bounds of the la statistical 
uncertainty as functions of the input WIMP mass. The dashed blue curves indicate the la band 
given with the parameters k and k' estimated analytically by Eqs. (35a) and (35b), whereas the 
solid red curves indicate the band given with the numerically estimated k and k'. Middle: dis- 
tributions of the reconstructed WIMP masses with analytically (dashed blue) and numerically 
(solid red) estimated k and k' . While the cyan (magenta) vertical lines indicate the median 
values of the simulated results (corresponding to the blue and red distribution histograms) and 
the horizontal thick (thin) bars show the 1 (2) a (68.27% (95.45%)) ranges of the results, the 
green vertical line indicates the true (input) WIMP mass of m x = 50 GeV. Right: reconstructed 
WIMP mass and the la statistical uncertainty bands given by the median values of the recon- 
structed Q VthIC with k and k! estimated analytically (dashed blue) and numerically (solid red), 
respectively. 28 Si and 76 Ge have been chosen as two target nuclei. 50 total events on average in 
each data set have been simulated. The input mass splitting has been set as 5 = 25 keV. Other 
parameters are as in Figs. 1. See the text for further details. 

the la band given with the parameters k and k' estimated analytically by Eqs. (35a) and (35b), 
whereas the solid red curves indicate the band given with the numerically estimated k and k' . 

It can be seen that, firstly, for input WIMP masses 20 GeV < m x < 70 GeV, while the 
analytically reconstructed WIMP mass (dashed blue) could be a bit overestimated, one could in 
principle reconstruct m x by the numerical method (solid red) pretty well. However, for heavier 
WIMP masses 70 GeV < m x < 500 GeV, m x reconstructed by both methods would be 
(strongly) underestimated. Nevertheless, the la upper bound could still cover the input (true) 
value and therefore offer at least an upper constraint on the mass of inelastic WIMPs. 

Meanwhile, in the middle frame of Figs. 8 we show the distributions of the reconstructed 
WIMP masses with analytically (dashed blue) and numerically (solid red) estimated k and k' . 
While the cyan (magenta) vertical lines indicate the median values of the simulated results 
(corresponding to the blue and red distribution histograms) and the horizontal thick (thin) bars 
show the 1 (2) a (68.27% (95.45%)) ranges of the results, the green vertical line indicates the true 
(input) WIMP mass of m x = 50 GeV. This plot shows that, for input WIMP masses 20 GeV 
< m x < 70 GeV, the reconstructed WIMP masses should be concentrated around the true 
(input) value with tails in the high-mass range. However, it has also been found that, once the 
WIMP mass is heavy {m x > 100 GeV), the distributions of the reconstructed m x could extend 
pretty widely (not only to the high-m x range, but also to the unphysical, negative range). 

On the other hand, in contrast to a wide spread of the distribution of the reconstructed m x , 
the distribution of the reconstructed Q VthTe has been found to in principle be concentrated and 
converged around the theoretical value. Therefore, in the right frame of Figs. 8 we show the 

reconstructed results, especially for heavier WIMP masses (m x > 40 GeV). 
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Figure 9: As in Figs. 8, except that the input mass splitting has been set as 6 = 10 keV (upper) 
and S = 40 keV (lower). Meanwhile, in the middle frames the input WIMP masses have been 
set as m x = 25 GeV (upper) and m x = 100 GeV (lower). 

reconstructed WIMP mass and the la statistical uncertainty bands given by the median values 
of the reconstructed Q VthIC with k and k' estimated analytically (dashed blue) and numerically 
(solid red), respectively. This plot shows that m x estimated by Eq. (21) with the median values 
of the reconstructed Q VthTC by using several data sets (with the same target nuclei) could indeed 
be (much) better than the median values of m x reconstructed by Eq. (21) with each single pair 
of data sets, especially for heavier input WIMP masses (m x > 100 GeV). 

As comparison, in Figs. 9 we consider the case of a smaller mass splitting of 5 = 10 keV 
(upper) and that of a larger one of 5 = 40 keV (lower). The left frame in the upper row of Figs. 9 
indicates that, while one could reconstruct the WIMP mass pretty well by both of the analytic 
and numerical methods for light WIMP masses (m x < 30 GeV), for heavier WIMP masses 
(m x > 50 GeV), the reconstructed m x would be (strongly) underestimated. Fortunately, a la 
upper bound as a constraint on the mass of inelastic WIMPs could still be given. 

The middle frame in the upper row of Figs. 9 shows that, (only) for available WIMP mass 
range (10 6 <5 < m x < a few x 10 6 <5), the distribution of the reconstructed m x could in principle be 
concentrated around the true (input) value, with however tails in the high- and negative-mass 
ranges. Moreover, the right frame shows that, for lighter mass splitting (e.g., 5 = 10 keV shown 
here), m x reconstructed by using the median values of Q Vthrc could also be (much) better than 
the median values of m x reconstructed with each single pair of data sets, up to even heavier 
WIMP mass of m x > 1 TeV. 

On the other hand, the lower row of Figs. 9 shows that, once the mass splitting 5 > 40 keV, 
except of the mass range between 40 GeV and ~ 120 GeV, the distribution of the reconstructed 



Figure 10: As in the upper frame of Figs. 9 (the input WIMP mass and the mass splitting has 
been set as m x = 25 GeV and 5 = 10 keV, respectively), except that 40 Ar and 136 Xe have been 
chosen as two target nuclei. 

WIMP mass could spread pretty widely and the reconstructed m x could anyway be (strongly) 
underestimated. 

Generally speaking, for the reconstruction of the WIMP mass by using the target combination 
of 28 Si and 76 Ge, our simulations show that, for mass splittings 5 < 40 keV, one could in prin- 
ciple reconstruct the WIMP mass in the range between 10 6 S and (3 — 5) x 10 6 <5 pretty precisely 
with a statistical uncertainty of ~ 30% (for m x ~ 10 6 <5) to a factor of ~ 2 (m x ~ 5 x 10 6 <5). 

Finally, in Figs. 10 we present our simulation results of the m x reconstruction by using the 
combination of 40 Ar and 136 Xe as our target nuclei. Here we show only the case of a light mass 
splitting of 5 = 10 keV. It has been found that, interestingly and a bit unexpectedly, by us- 
ing the target combination of 40 Ar and 136 Xe the strong underestimation of the reconstructed 
WIMP mass in the high-mass range could be (strongly) alleviated with a much smaller statis- 
tical uncertainty (left frame) and the distribution of the reconstructed m x would also be more 
concentrated (middle frame). 

3.4 Determining d 

In this subsection, we present the simulation results of the reconstruction of the second key 
properties of inelastic WIMPs, i.e., the mass splitting 5. 

In the left frame of Figs. 11, we show the reconstructed mass splitting estimated by Eq. (22) 
and the lower and upper bounds of the la statistical uncertainty estimated by Eq. (24) as 
functions of the input mass splitting for the case of m x = 100 GeV by using 28 Si and 76 Ge as two 
target nuclei with 50 total events on average in each data set. The dashed blue curves indicate 
the la band given with the parameters k and k' estimated analytically by Eqs. (35a) and (35b), 
whereas the solid red curves indicate the band given with the numerically estimated k and k'. 

It can be seen that, firstly, due to the maximal kinematic cut-off energy Q max ,kin in Eq. (42b), 
one could reconstruct the mass splitting only up to 5 ~ 50 keV (for the WIMP mass of ~ 100 
GeV). However, the middle frame of Figs. 11 shows that, once this reconstruction is achievable, 
the mass splitting could in principle be reconstructed pretty well: not only the median values of 
6 reconstructed both analytically and numerically could match the true (input) one, but also the 
distribution of the reconstructed 5 would be pretty concentrated. Moreover, the middle frame 
of Figs. 11 shows also that a 3a to 5a (> 99%) CL of the positivity of the reconstructed mass 
splitting (as the confirmation of the inelastic WIMP scattering) could in principle be identified. 

Finally, as in Figs. 8, in the right frame of Figs. 1 1 we show the reconstructed mass splitting 
and the la statistical uncertainty bands given by the median values of the reconstructed Q Vthie 



Figure 11: Left: reconstructed mass splitting and the lower and upper bounds of the la statistical 
uncertainty as functions of the input mass splitting. The dashed blue curves indicate the la 
band given with the parameters k and k' estimated analytically by Eqs. (35a) and (35b), whereas 
the solid red curves indicate the band given with the numerically estimated k and k'. Middle: 
distributions of the reconstructed mass splittings with analytically (dashed blue) and numerically 
(solid red) estimated k and k' . While the cyan (magenta) vertical lines indicate the median 
values of the simulated results (corresponding to the blue and red distribution histograms) and 
the horizontal thick (thin) bars show the 1 (2) a (68.27% (95.45%)) ranges of the results, the green 
vertical line indicates the true (input) mass splitting of 5 = 25 keV. Right: reconstructed mass 
splitting and the la statistical uncertainty bands given by the median values of the reconstructed 
Q VthTe with k and k! estimated analytically (dashed blue) and numerically (solid red), respectively. 
28 Si and 76 Ge have been chosen as two target nuclei. 50 total events on average in each data set 
have been simulated. The input WIMP mass has been set as m x = 100 GeV. Other parameters 
are as in Figs. 1. See the text for further details. 

with k and k! estimated analytically (dashed blue) and numerically (solid red), respectively. Not 
surprisingly, 5 estimated by Eq. (22) with the median values of the reconstructed Q Vthrc could be 
so good as the median values of 5 reconstructed by Eq. (22) with each single pair of data sets; 
the former could have little bit smaller statistical uncertainties. 

As comparison, in Figs. 12 we consider the case of a smaller WIMP mass of m x = 50 GeV 
(upper) and that of a larger one of m x = 250 GeV (lower). In contrast to the reconstruction of 
the WIMP mass shown in Figs. 8 and 9, our simulations show here clearly that, for a rather 
small mass splitting 5 ~ a few (tens) keV, one could always reconstruct 5 pretty well, with only 
a bit larger statistical uncertainty for large WIMP masses. More importantly, a clear 3a to 5a 
(> 99%) CL of the positivity of the reconstructed mass splitting could always be identified. 

In fact, it has been found that, up to a WIMP mass of ~ 1 TeV, one could in principle 
always reconstruct the mass splitting pretty precisely with a deviation of < 20% (analytically) 
or even < 10% (numerically) and a statistical uncertainty of ~ 15% to ~ 50% (analytically) or 
~ 10% to a factor of ~ 2 (numerically). Note that, as shown in Figs. 11 and 12, the statistical 
uncertainty on the analytically reconstructed mass splitting increases with increasing the mass 
splitting, whereas that of the numerically reconstructed one is approximately the same for all 
values of the mass splitting. 

3.5 (5 = (elastic WIMP— nucleus scattering) case 

In this subsection, we consider the special case of the zero mass splitting 5 = 0, i.e., elastic 
WIMP-nucleus scattering, in order to demonstrate the ability of distinguishing the inelastic 
WIMP scattering scenarios from the elastic one. 
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Figure 12: As in Figs. 11, except that the input WIMP mass has been set as m x = 50 GeV 
(upper) and m x = 250 GeV (lower). Meanwhile, in two middle frames the input mass splittings 
have been set the same as 5 = 25 keV. 

In the first column of Figs. 13, we show the theoretical recoil spectra (dashed-double-dotted 
cyan), the measured energy spectra (dotted magenta histograms) as well as the analytically 
and numerically reconstructed energy spectra with four different target nuclei. It can be seen 
that, for this special case with all four simulated detector materials, while the analytically 
reconstructed spectra are still somehow peaky (with non-zero k' as well as non-zero Q Vthie ), the 
numerical iterative process could now indeed offer a better (and better) reconstruction of the 
recoil spectrum with much smaller k' and Q Vthia (see the third and forth columns). Meanwhile, the 
distributions of the reconstructed k' and Q VthIC indicate consistently that, although the median 
values of the (numerically) reconstructed k' and Q Vthia would be non-zero, the most frequently 
observable values of k' and Q VthIC could however be negligibly small! 

Moreover, it has also been found that, importantly, for this special S = (elastic scattering) 
case, although the reconstructed mass splitting would be a little bit non-zero (positive), one 
would observe simultaneously non-physically a negative reconstructed WIMP mass (m x < 0); 
the larger the true (input) WIMP mass, the larger the absolute value of the reconstructed one. 
This unique observation could in turn help us to confirm or rule out the inelastic WIMP-nucleus 
scattering scenarios down to a (very) small mass splitting (5 < a few keV). 

Here we would like to emphasize that, even though the measured recoil energy spectra (his- 
tograms shown in the first columns of Figs. 3 and 13) look almost the same (and thus would 
be difficult to be distinguished from each other by conventional data analysis), our model- 
independent (analytic and numerical) reconstructions of the recoil spectrum as well as the esti- 
mations of Q VtbTC could give clearly different results (cf. the third and forth columns of Figs. 3 
and 13). 




( 136 Xe) 

Figure 13: As in Figs. 1 and 2, except that the zero mass splitting 5 = has been used. 

4 Summary and conclusions 

In this paper we study direct Dark Matter detection experiments for the inelastic WIMP-nucleus 
scattering framework and develop model-independent methods for not only reconstructing the 
measured recoil energy spectrum, which can in turn be used for identifying the inelastic WIMP 
scenarios, but also determining the WIMP mass and the mass splitting. 

At the beginning of this work, we followed our earlier work [27] to derive the formula for 
reconstructing the one-dimensional velocity distribution function of inelastic WIMPs. However, 
since for inelastic WIMPs, there are two unknown parameters: the (degenerate) mass m x and 



the tiny mass splitting 5, not only our formula for reconstructing fi(v), but also the method 
for determining the WIMP mass introduced in Ref. [28], can not be used directly. Hence, we 
turned to develop a new procedure for determining the WIMP mass as well as the mass splitting 
model-independently. 

For this aim, we introduced a two-parameter exponential ansatz for reconstructing the mea- 
sured recoil energy spectrum as well as determining the characteristic energy corresponding to 
the threshold (minimal required) velocity of incident inelastic WIMPs, which can produce recoil 
energy at all. In this process, not only the additional fitting parameter for reconstructing the 
recoil energy spectrum is approximately proportional to the squared mass splitting, but also the 
characteristic energy is directly proportional to the mass splitting (for a fixed, non-zero WIMP 
mass). Thus these two quantities could be good indicators for identifying inelastic WIMP- 
nucleus scattering scenarios. Our numerical simulations show that, not only for large (~ a few 
tens keV) but also for (very) small (~ a few keV) mass splittings, our model-independent re- 
construction could in principle indeed identify the positivities of the additional fitting parameter 
and the characteristic energy of the recoil energy spectrum with a 3a to 5cr confidence level up 
to a WIMP mass of ~ 1 TeV and a mass splitting of ~ 200 keV. 

Meanwhile, for analytically reconstructing the recoil energy spectrum, one has to assume that 
the minimal cut-off energy could be negligible and the maximal one is large enough. Although 
these are achievable in most currently running and next-generation experiments, the kinematical 
minimal and maximal cut-off energies would limit the available ranges of the mass (splitting), 
especially for the use of light target nuclei. Therefore, we introduced in this work a numeri- 
cal iterative procedure for correcting the analytic estimations of the fitting parameters of the 
measure recoil energy spectrum. Our simulations show that, basically this numerical correction 
could indeed offer the better reconstructions with smaller statistical uncertainties. However, the 
statistical fluctuation due to the use of only a few tens events causes the divergency problem of 
the distribution of the characteristic energy. For larger mass splittings (5 > 40 keV), the ana- 
lytic reconstructions seem to be more reliable; the distributions of the reconstructed spectrum 
fitting parameters and the characteristic energy offered by the numerical iterative procedure 
could have a (much) longer tails in high- (and even low-)value ranges. This problem would be 
worse for the use of light target nuclei; results offered with heavy target nuclei could still be 
used as auxiliary. 

Moreover, we considered several different target nuclei. It has been found that, once the mass 
splitting is small (5 < 40 keV), all materials could identify inelastic WIMP scattering pretty 
well. In practice, light target nuclei, e.g., 28 Si and 40 Ar, could even work better, due to relatively 
higher values of the second fitting parameter and the characteristic energy of the reconstructed 
recoil spectrum. However, for larger mass splittings (5 > 50 keV), limited by the kinematical 
minimal and maximal cut-off energies due to the Galactic escape velocity of halo WIMPs, one 
would need heavy target nuclei, e.g., 76 Ge and 136 Xe, for the identification of inelastic WIMPs. 

Furthermore, for the reconstruction of the WIMP mass with the target combination of 28 Si 
and 76 Ge, our simulations show that, for mass splittings 5 < 40 keV, one could in principle 
reconstruct the WIMP mass in the range between 10 6 5 and (3 — 5) x 10 6 <5 pretty precisely with a 
statistical uncertainty of ~ 30% (for m x ~ 10 6 5) to a factor of ~ 2 (m x ~5x 10 6 5). Meanwhile, 
we found also that the WIMP mass estimated with the median values of the reconstructed 
characteristic energy by using several data sets (with the same target nuclei) could indeed be 
(much) better than the median values of the WIMP mass reconstructed with each single pair of 
data sets, especially for heavier input WIMP masses (m x > 100 GeV). Moreover, it has been 
found that, once the mass splitting is pretty light (5 ~ 10 keV), by using the target combination 
of 40 Ar and 136 Xe, the strongly underestimated WIMP mass in the high-mass range with the 



Si-Ge combination could be (strongly) alleviated with a much smaller statistical uncertainty; 
the distribution of the reconstructed m x would also be more concentrated. 

On the other hand, due to the limitation by the maximal kinematic cut-off energy and 
the required use of two (one heavy and one light) target nuclei, one could only reconstruct 
the mass splitting less than a few tens keV. Nevertheless, in the reconstructable range, one 
could in principle always reconstruct the mass splitting pretty precisely with a deviation of 
< 20% (analytically) or even < 10% (numerically) and a statistical uncertainty of ~ 15% 
to ~ 50% (analytically) or ~ 10% to a factor of ~ 2 (numerically) up to a WIMP mass of 
~ 1 TeV. Whether and how to use this (pretty) precisely reconstructed mass splitting as a priorly 
determined parameter for further reconstructions of the mass as well as the one-dimensional 
velocity distribution of inelastic WIMPs will be investigated in the future. 

Finally, we consider the special case of the zero mass splitting (elastic WIMP-nucleus scat- 
tering). It has been found that, firstly, for this special case the numerical iterative procedure 
could offer an estimation of (almost) zero of the characteristic energy of the recoil spectrum. 
In addition, the detailed analysis of the distributions of the reconstructed results shows that, 
although the median values of the (numerically) reconstructed characteristic energies would be 
a little bit non-zero, the most frequently observable values could however be negligibly small! 
Secondly, it has also been found that, although the reconstructed characteristic energy as well 
as the reconstructed mass splitting would be a little bit non-zero (positive), one would observe 
simultaneously non-physically a negative reconstructed WIMP mass; the larger the true (input) 
WIMP mass, the larger the absolute value of the reconstructed one. This unique observation 
could in turn help us to confirm or rule out the inelastic WIMP-nucleus scattering scenarios 
down to a (very) small mass splitting (5 < a few keV). Moreover, we would like to emphasize 
that our model-independent (analytic and numerical) reconstructions of the recoil spectrum as 
well as the estimations of Q VthIC could give clearly different results between the inelastic scattering 
case of a small, but non-zero mass splitting and the elastic one. 

In summary, as complementarity and extension of our earlier work on the development of 
(model-independent) methods for reconstructing properties of Galactic WIMPs by using direct 
DM detection data, we introduce in this work new model-independent approaches for identifying 
inelastic WIMP-nucleus scattering as well as for reconstructing the mass and the mass splitting of 
inelastic WIMPs simultaneously and separately. Our results show that, with a few tens observed 
WIMP signals (from one experiment), one could already distinguish the inelastic WIMP scenarios 
from the elastic one. By using two or more data sets with positive signals, the WIMP mass and 
the mass splitting could even be reconstructed. Hopefully, these could help our experimental 
colleagues to not only distinguish different frameworks of DM/particle physics, but also further 
constrain the parameter space in (various) extensions of the Standard Model of particle physics. 
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A Expressions for the two— parameter exponential ansatz 

In this section we give detailed derivations and expressions needed for numerical estimations of 
different moments of the two-parameter exponential ansatz. 

A.l Estimator of r given in Eq. (46) 

We start with the inequality of arithmetic and geometric means: 

\(ax + -\ > Jax-- = y/aVb. (Al) 

Then we can define that 

ax + - = 2y/aVb cosht > 2yW&, Vie [-00,00]. (A2) 

x 

Differentiating the both sides, one can get 

(a — ^- ) alx = — ( ax — — ) dx = 2\fa\fb sinht alt . (A3) 
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On the other hand, from the definition (A2), we have 
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i.e., 



2y/aVb sinh t = ax — — . (A6) 
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Comparing the above equation with Eq. (A3), we have 
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Therefore, we can obtain that 
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Here, firstly, from Eq. (A2), one can find that, once x = yb/a, cosht = 1, and thus t = 0. 
Secondly, for the last line we have used the integral formula for the modified Bessel function of 
the second kind: 
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Finally, by differentiating with respect to a, we can obtain an analytic form for the denominator 
of ro in Eq. (31) for the case of a negligible minimal cut-off energy and a (very) large maximal 
one as 
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A. 2 Moments of the two— parameter exponential ansatz 

Since we have found that 
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by differentiating with respect to a, one can get 
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Similarly, by differentiating with respect to 6, we have 
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Therefore, by differentiating with respect to a, one can get 
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Similarly, with respect to b, we have 
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On the other hand, by setting x = 1/y, one can find that 
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for A = 1/2, 1, 3/2, 2, • • •. Hence, one can obtain Eqs. (Al3b), (Al3c), (Al8b) and (A18c) 
from Eqs. (A12) and (Al3a) and Eqs. (A17) and (Al8a) by exchanging a <-> b and x ■<-)■ 1/x and 
using 

erf (— x) = — erf(x) . 



B Solving the fitting parameters k and k' numerically 

In this section, we describe our numerical iterative procedure for solving the fitting parame- 
ters of the measured recoil energy spectrum based on the analytically estimation of these two 
parameters. 




Figure Al: Sketch of the numerical procedure for giving a linear equation of k and k! by using 
the minus-one-half moment of the two-parameter exponential ansatz. 

We start with the point (fcanaj^Lia)) which has been estimated by Eqs. (35a) and (35b). 
Check at first whether (l/y/x) (k ana , k' ana ; x) - (l/VxY = F_ 1/2 (k = k aDa ,k' = k' ana ) > 0. 



Here (l/y/xj (k, k'; x) is the function of k and k! given by Eq. (A17) and 
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can be estimated from measured recoil energies directly. Note that the integral on the right-hand 
side of Eq. (A21) has an upper (lower) bound of Q max (Q m in) and can be estimated numerically 
by setting, e.g., fcW* = fc^. 

Without losing the generality, we assume that F_ 1 / 2 (k !ma , k' aDa ) > 0. Then set (k ana , k' ana ) as 

(k mui , k' min ) . Note that the function (\/y/x) (k, k'; x) Q should decrease monotonically as k 

^ ' Qmin 

and/or k' increases. As sketched in Fig. Al, by, e.g., fixing k' = k' min and then checking whether 
F_i/ 2 {k, k') > by setting k = k min + n (/c min /10), n — 1, 2, 3, • • •, one should be able to 
find two auxiliary points (k max , k f min ) and (k min , k' nmx ) at which the function F_ 1 / 2 {k,k') has a 
different sign (< corresponding to our current assumption) Now solve A; sol G [k min , k max ] and 
*4>i £ [^mini^max]! which satisfy that F_ 1/2 (k soh k' min ) < 1(T 5 (l/^/x)* and F_ 1/2 (k min , k' sol ) < 

io- 5 (i/^r. 

Since the intersection boundary of F_i/ 2 (k, k') = on the k — k' plane is almost a straight 
line and the analytic estimates (k ana , k ana ) should in principle be pretty close to the numerical 
solution (A; num , k' nnm ), we approximate the equation of the boundary to a linear equation by 
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Then the numerical solution of (k, k') can be given by 
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Here we define 
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Here k^l m 1 indicate the numerical estimates of the fitting parameters in the first round. One 
could process the whole numerical procedure iteratively. Remind however that statistical fluc- 
tuation could cause a divergency problem and the results from the later rounds might be worse 
than that from the first or second round. 
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C Expressions of the derivatives of Q Vthre , k and h! 

Firstly, differentiating the expression (38) for solving the characteristic energy Q Vthie with respect 
to k, one find that 
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Similarly, differentiating the expression (38) with respect to k', one has 
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Thus, it can be found that 
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Remind that, as Eqs. (36) and (38), Eqs. (A29a) and (A29b) can be used for both of the 
analytically and numerically estimated k and k' . 



C.l For the analytic estimates 

For the analytic estimate of k given by Eq. (35a), we have 
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Similarly, from the expression (35b) of k', we have 
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C.2 For the numerical solutions 



Firstly, from Eq. (43), we have 
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Therefore, by using definitions (A17), (Al8a) and (Al8b), one can get that 
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